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04 , We propose a method to calculate dispersion interactions in a system composed of a one dimen- 

sional layering of finite thickness anisotropic and optically active slabs. The result is expressed 
' within the algebra of 4 x 4 matrices and is demonstrated to be equivalent to the known limits of 

, isotropic, nonretarded and uniaxial dispersion interactions. The method is also capable of handling 

■ dielectric media with smoothly varying anisotropy axes. 
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. der Waals or dispersion torques that have the same origin as the more common and ubiquitous long range forces. 
' Following this seminal contribution the effect of dielectric anisotropy in the dispersion interactions between two 
semiinfinite birefringent media has been studied theoretically in several contexts 0, Q and has been given the most 
thorough analysis by Barash and coworkers (1, The outcome of these theoretical advances was that apart from 
the dispersion interactions that are present for any dielectric media, two birefringent slabs will also exhibit dispersion 
torques acting to align them. The magnitude of this torque is of course small but should nevertheless be detectable. 
Inspite of these theoretical advances there are as yet no experimental verifications of the theoretical predictions though 
there is plenty of activity in this research field [7[ ■ It goes without saying that the dispersion torques should play a 
Q ' very important role in the context of micro and nano-mechanical actuators that could in principle transduce not only 
O . rectilinear motion but also rotation, or even couple and convert the two. 

' In the present work we will formulate the theory of dispersion interactions in a system composed of a one dimen- 
I \ sional layering of finite thickness anisotropic dielectric slabs. In this way we generalize the already existing theory 
^ ■ of dispersion interactions in a multilayer system of isotropic dielectric materials [8]. In this approach the secular 
' determinant of the modes, which enters the free energy of the dispersion interactions, can be obtained as one of the 
. elements of the transfer matrix of the EM field modes propagating through the system. The complete transfer matrix 
of the complete stratified system can be decomposed into a product of two separate matrices, a diagonal propagator 
, matrix and a symmetric discontinuity matrix. This method of matrix decomposition is extremely elegant and time 
' saving for complicated one-dimensional stratifications and has been since used profitably in a variety of contexts . 
I The matrix approach introduced in the cited work has been formulated only for the case dispersion interactions 
between dielectrically isotropic media. For a complete solution of the anisotropic case one would need to extend it to 
solutions of a complete set of Maxwell equations with tensorial response functions. Indeed, we will show below that 
it can be extended and generalized to solve the dispersion interactions problem across stratified non-isotropic media. 
While the isotropic case can be formulated with the help of an algebra of 2 x 2 matrices, the general anisotropic case 
^ . can only be implemented within an algebra of 4 x 4 matrices and is thus notoriously more difficult to handle, but 
d ' still simpler than would follow from alternative approaches [1]. Nevertheless we are able to extract some informative 
limiting cases that show the power and book-keeping elegance of this approach. 
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II. FORMALISM 



The propagation of electromagnetic waves in layered anisotropic media was described already by Berreman who 
introduced a convenient and efficient formalism that allows for a transparent formulation of a notoriously extremely 
complicated problem In what follows we basically follow this approach. If z is the direction perpendicular to the 
translational plane of symmetry of the problem, so that the dielectric slabs of anisotropic material are layered in this 
direction, the propagation of the single frequency electromagnetic field in such a system can be Fourier-decomposed in 
the X — y plane perpendicular to the z direction with wavenumber vector Q and described using the four dimensional 
vector of the field components 



^(z;g) 



E,{z;Q) 
Hy{z;Q) 
Ey{z;Q) 
-H,{z;Q) 



(1) 



The Fourier field components evolve as a function of the coordinate z via a homogeneous system of linear differential 
equations that stem directly from the properly formulated Maxwell equations 



9 , > 

7r^ = — 

oz c 



(2) 



where D = T){e{uj, z), ii{uj, z), p{uj, z), p'{uj, z),Q,ll!) is a 4 x 4 complex matrix that depends on the frequency and z 
coordinate dependent 3x3 tensors of the electric and magnetic permeabilities e and fi, the 3x3 tensors of optical 

activity p, p', the twodimensional vector Q in the x — y plane that gives the wavenumber in the plane of symmetry 
and the chosen angular frequency oj. 

The details on how to construct D are given in [ll|, with the main results repeated for consistency here. A 6 x 6 
matrix M is constructed from the material response tensors as 



M = 



e P 
p' ^ 



(3) 



An auxiliary matrix a is then calculated as 



with 



03,1 = 


[^6,1^3,6 - M3,lM6,6]/rf, 


(4) 


03,2 = 


[(M6,2 - CQ/C^)M3,6 - M3,2M6,6]/d, 


(5) 


03,4 = 


[M6,4M3,6 - M3,4M6,6]/d, 


(6) 


03,5 = 


[M6,5M3,6 - (M3,5 + cQ/t^)Af6,6]/rf, 


(7) 


06,1 = 


[M6,3M3,i - M3,3M6,i]/d, 


(8) 


06,2 = 


[M6,3M3,2 - M3,3(M6,2 - cQ/w)]/d, 


(9) 


06,4 — 


[M6,3M3,4 - M3,3M6,4]/rf, 


(10) 


06,5 = 


[M6,3(M3,5 + cQ/t^) - M3,3Af6,5]/rf, 


(11) 




d = M3,3M6,6 - Af3,6M6,3- 


(12) 



Without any loss of generality the wavenumber Q is assumed to lie in the x direction and we have to rotate the 
coordinate frame accordingly. The elements of a 4 x 4 matrix S are then defined as 



^1.1 


= Afi,i- 


1- A'/i,3a34 - 


V Mi^6a6,i, 


(13) 


<5'l,2 


= Mi,2- 


h Mi^3a3^2 - 


1- Mi,6a6,2, 


(14) 


'5'l,3 


= Ml 4 - 


V Mi^3a3^4 - 


V Mi_6a6,4, 


(15) 


51,4 


= Mi,5- 


1- Mi^3a3^5 - 


1- ^1^606,5, 


(16) 


<S'2,1 


= M2,l - 


V M2,3a34 - 


V {M2fi - cQ/uj)ae^i, 


(17) 


<5'2,2 


= M2,2- 


V Mi,za.za - 


\- (Af2,6 - cQ/uj)ae^2, 


(18) 


<5'2,3 


= M2,4- 




\- (Af2,6 - cQ/uj)ae^4, 


(19) 
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5*2,4 = ^2,5 + M2,3a3,5 + (M2,6 - cQ/w)a6,5, 

5*3,1 = M4,i +M4,3a3,l + M4,6a6,l, 

53.2 = ^^4^2 + -^"4,303, 2 + ^4^606, 2, 

53.3 = M44 + M4^3a3,4 + Af4,6a6,4, 

53.4 = M4^5 + M4^3a3,5 + M4,6a6,5, 

54.1 = Af54 + (M5,3 + c(5/Lj)a3,i + Ms^goe,!, 

54.2 = Af5^2 + (A^5,3 + c(3/a;)a3,2 + -^5,606,2, 

54.3 = M<SA + (^5,3 + cQ/uj)a3^4 + M5,6a6,4, 

54.4 = Af5,5 + (Af5,3 + cQ/uj)a3^5 + M5,6a6,5, 

so that the sought for matrix D is then given as 

D = H J S H^^ 

with a matrix corresponding to a permutation of components 



H 



10 
1 
10 
10 



(20) 
(21) 
(22) 
(23) 
(24) 
(25) 
(26) 
(27) 
(28) 



(29) 



(30) 



and a transformation of the form 








1 

-1 



0-100 
10 



(31) 



If the matrix D is constant within an interval [z, z + h], then the propagation of EM field modes within that interval 
can be written in terms of a matrix exponent 



^p{z + h) ^P{h)^{z), P(/i)=exp 



-D h 



(32) 



Here the matrix P is a constant within that interval but does depend on the length of the interval, h. A matrix 
exponent of the form above can be calculated by performing a diagonalisation of D 



P{h) = RK{h) K-\ Kjj ^ exp (^—dj hj 



(33) 



where K is a diagonal matrix, dj are the eigenvalues of the matrix D and the rotation matrices R and correspond 
to the left and right eigenvectors of D. The matrix performs a mapping from the vector i/j of the electromagnetic 
field components into the vector of two left and two right propagating plane wave components that we denote by 77. 

Propagation of EM waves across a stratified set of — 1 layers of finite width, that are bounded by semi-infinite 
half-spaces from both directions, can be described by taking a product 



Pjv(/ijv) Pjv-i(/ijv-i) . . . Pi(/ii) Poiho) 



(34) 



or 

P RArKAr(/l„)R^"^ RjV_ iKjV- 1 (/in- 1 )RjV- 1 ■ ■ • R'oKo(/io)Ro ^- (35) 

where the indices and N correspond to the left and right half-spaces. This is indeed very similar to the case treated 
in but in that case the dielectric media are assumed to be isotropic. 

A propagator of the form Eq. [35] maps the fields from a point that is at some arbitrary depth ho in the left half-space 
to a point that is at another arbitrary depth /ijv in the right half-space. The propagators Pq and Pat, however, do 
nothing more than rotate the phases of the incoming and outgoing plane waves. 

Furthermore, the propagator P maps the components of the electric and magnetic fields, whereas it is the mapping 
of plane wave amplitudes that is of interest for dispersion interactions. A mapping of plane wave amplitudes from 
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one side of the stratified system of slabs to the other is obtained by the truncation of P on both the left and the right 
sides of the expression to obtain 



P — R-at"^ riAf-lKjv_i(/lAr_i)RjY"'_i ■ • • RlKi(/li)Rj ^ Ro 



(36) 



In order to propagate the vector of plane wave amplitudes from the external boundary of the leftmost layer to the 
boundary of the rightmost one than has 



(37) 



This propagator does not depend on the arbitrary thicknesses Hq and hj^ anymore. The only information that the 
bounding half-spaces contribute to the propagator is their respective transformations from the electromagnetic fields 
into plane waves as given by the matrices Rq, R^v^- 

With these expressions it is now possible to tackle the problem of dispersion interactions in such media in analogy 
with j9|. The dispersion interaction free energy F per area A is given as 



F/A 



CO ' „ 

n=0 •' 



(2^)2 



lnC(iC„). 



(38) 



The prime in the summation over the Matsubara thermal frequencies denotes that the n = term is to be taken with 
the weight of 1/2. Here C{uj) is an expression whose zeros give the eigenfrequencies of bound states. It is evaluated 
as a function of the imaginary Matsubara frequencies ^„ = 2TTnkT/h. In order to evaluate the above interaction free 
energy we first need to calculate C(a;) in terms of the propagator equation Eq. [ST] 

We define a state as bound if it is exponentially decreasing to both —oo and to +00. Let the four components of 
the plane wave amplitude vector rj for each layer be sorted with respect to their real parts such that the exponentially 
increasing solutions have indices 1, 2 and the decreasing ones have indices 3, 4. In order for the state to be bound, the 
exponentially decreasing components of rjQ and the exponentially increasing components of rjN in the corresponding 
half-space must be equal to 0. This gives for the wave propagation equation [37] corresponding to bound states the 
form 



,f ) = 



Pi Pn 
Pin Piv 



(0) 

Vi 

(0) 







(39) 



(40) 



that can be clearly reduced to a 2 x 2 homogeneous system of equations 

This system of equations allows us to identify the bound state condition that gives their eigenfrequencies as 

C = detP/ = 0. (41) 

When evaluating C at imaginary frequencies in equation (|38p , the determinant det P/ will in general have a complex 
value. We may take its absolute value as the appropriate contribution to the free energy, since the values of C for 
the opposite wavenumbers Q and — Q are complex conjugate and therefore the imaginary values of their logarithms 
cancel. 

As the free energy expression is indeterminate up to an additive constant, this means that the bound state condition 
is undetermined up to a constant scaling factor. It is natural to choose such a scaling that further splitting of the 
half-spaces into sublayers with equal properties does not alter the results. This will also ensure that the dispersion in- 
teractions only depend on the spatial variations of the electromagnetic response tensors and not on their homogeneous 
values. This is acheived by shifting the eigenvalues of the D matrix by 



di — di ~ d, d = {di + d'l)!'! 



(42) 



and use these when evaluating the propagators K^. 

There is still an additional freedom of scaling factors due to the eigenvector normalisations. Since the matrix D 
is in general not hermitian, its left and right eigenvectors are not complex conjugate and there is no unique way to 
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normalise them individually. Care must therefore be taken when comparing systems whose terminating half spaces 
are being rotated with respect to each other. A common reference point (such as the case of large separation) must 
be calculated in order to be able to determine the free energy shift due to eigenvector rotations. 

One particular reference point for any system can be a corresponding system in which all the layers between the 
two limiting half-spaces are substituted by vacuum, and the limit of the distance between the two half-spaces (the 
vacuum distance) is taken to infinity. This configuration presents no interaction between the two half-spaces, and 
therefore the calculated free energy contribution of this configuration can be considered as the sought shift. For this 
corresponding system, the matrix Pc is given by 



P 



Rn 



where Poo corresponds to the vacuum propagator for the electromagnetic waves in the limit of long distances, 
culating it using the expression and the shift ([i^ gives 



(43) 
Cal- 



1 

\ 









1 

2 

_£ 
2 






_ J_ 

2p 

I 
2 



(44) 



The free energy contribution for this corresponding system is then calculated via pd|) and and (|38p . In this way one 
can determine the free energy shift due to eigenvector rotations and properly normalize the dispersion interaction free 
energy. This concludes our derivation of the general formalism. In what follows we will consider a few interesting 
limiting cases. 



III. ISOTROPIC CASE 

When the response tensors are isotropic, the general approach introduced above should reduce to the formalism 
introduced in Let us assume e — e/, /i = /x/ and p = p' = 0. The matrix D for this case reduces to the simple 
form 



D 




e 





(45) 



where ^ = — icj is the imaginary frequency. The eigenvalues of such a matrix are doubly degenerate and come in pairs 
of c?i = ±i^, where 



c2Q2 



(46) 



such that D = R diag {di} R ^ can be written in the form 



D 



-V V Q 

e e 

-/I 

Q Q V V 

















V 














— V 














V 



2v 2e 





2u 2v 
1 1 



0^^ 



(47) 



The matrices in the above expression are block diagonal, so the four dimensional problem decomposes into a pair of 
two dimensional ones that correspond to the TE and TM modes of the electromagnetic propagation. Let us elaborate 
this further. 

After shifting the eigenvalues in order to remove the exponential divergence, the propagator P is then given as 



R 



10 

cxp(-2p) 

10 

exp(-2p) 



R 



(48) 
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where p = in order to conform with the notation in Q. The transition of eigenmodes between two subsequent 
layers is given by the matrix 



£i + l 



Pi 



Pi+1 
Pi 



pi+1 Ci + l Pi+1 

















Mi _j_ Pi Pi I 

Pi+1 Pi+1 

Pi I pi Pi _|_ 

Pi+1 Pi+1 Pi+1 



y _P^ 

p.+l 
Pi 

p.+l J 



(49) 



If we further spht the above matrix into separate TE and TM diagonal blocks, and apply two arbitrary scaling factors 
due to the eigenvector normalisations by dividing the diagonal matrices with their diagonal terms, we obtain 



with 



for the TM and 



for the TE block. By also denoting 



1 -Ai 

-A,; 1 



Pi+lfJ-i — PifJ-i+1 
Pi+lP-i + PifJ'i+1 



1 

exp{-2pi) 



we can obtain the free energy contributions as the {1,1} matrix element of a matrix product, 

Cte, TM — [GnTn ■ ■ ■ GiTiGo]^ii-^ , 



(50) 



(51) 



(52) 



(53) 



(54) 



which agrees exactly with the result derived in Q . 

In the full 4x4 formulation there does not exist a single common scaling factor for the matrix R which would 
exactly map it to the 2x2 formalism of [9]. Scaling factors are, however, only related to the shift in the free energy, 
which is itself undetermined up to a constant. 



IV. THE SINGLE SLAB ANISOTROPIC CASE 



The nonretared case of two semi-infinite anisotropic dielectric media separated by a slab of a third anisotropic 
medium, where for all the media one of the anisotropy axes coincides with the z axis, was derived in [1]. This result 
can be summarized as follows. 

Let the dielectric tensors for the first half space, the slab and the second half space be given by e^, ^-^d £3, 
respectively. Each tensor can be written in the form 



(55) 



e-i,i 






















ei,3 



Each of the materials is rotated around the z axis such that 



where 



O, = 



cos(6',) -sin(6li) 
sin(6'i) cos(6'i) 
1 



(56) 



(57) 
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From a given deilectric tensor the quantity 

sin^(a) (58) 

£1,3 £1,3 

may be calculated. From this 

/3fe,^.,^)=QG(e„0,-v?) (59) 

is defined, where ip is the direction of the plane parallel wavevector for which the free energy contribution is being 
evaluated. Further 

a = '^^^^J^^^'l^'^l (60) 



and 



are introduced, from which 



e2,3^(e2:^2,<^) 

e3,3/3(e3,6>3,y)) 
£2,3/3 (£2 > ^2, <^) 



(61) 



= (62) 
(a + l)(6+l) 

is calculated. The free energy contribution for a fixed size of the plane parallel wave vector Q and its angle tp is given 
by 

Cp(ei, £2, £3^ ^1. ^2, 03, Q, V', i) = 1 - A2 exp {-2QGie^, 62 - p)L) , (63) 
and the total free energy then follows as 

F/A = -^Y. I '^^ I InCp Q dQ. (64) 

ji=0 



"'0 



While we were not able to demonstrate analytical equivalence of our approach with the above result, numerical 
computations nevertheless agree perfectly , up to a constant shift. As already explained, the constant free energy shift 
in our formulation is a result of arbitrary eigenvector normalization, and the results must be compared to some chosen 
point, which in this case is the infinite separation of the two materials. Note that this constant energy shift makes no 
contribution to the interaction forces, given by the derivative of the free energy. 

For comparison purposes, we choose the free energy contributions from our method, C, in the limit of inifinite 
c, and the corresponding nonretarded expression, Cp, calculated for the dielectric response tensors at some given 
imaginary frequency. These are arbitrarily chosen to be ei^i = 5, ei^2 = 10, £1,3 = 15, £2,1 = 1, £2,2 = 3, £2,3 = 2, 
£3,1 ~ 18, £3,2 — 12, £3,3 = 6. The second medium is rotated by the angle 62 = 2 with respect to the first one and the 
third one by ^3 = 1 with respect to the first one. The angle for the wavevector was chosen as tp = 3. The free energy 
contributions for both calculations as a function of distance are given in Fig. [T] For our method, the large distance 
free energy result was subtracted from all the results in order to eliminate the shift in the free energy contribution. 
The match is exact for this set of parameters as well as for any other one might choose. This strongly supports the 
general conclusion that the non-retarded form of our calculation coincides exactly with the calculation derived in [l[ . 

A calculation for the fully retarded dispersion interaction of the same configuration is given in , with the limitation 
that the materials are uniaxial and that the axis of anisotropy lies in the symmetry plane of the problem, with the 
medium between the two slabs being isotropic. The result is convoluted and lengthy, and will thus not be repeated 
here (for details see [1, 0, @])- Using the same notation as before, the arbitrarily chosen values are £i_i = 5, £1.2 = 10, 
£1,3 = 10, £2,1 = 2, £2,2 = 2, £2,3 = 2, £31 = 6, £3,2 = 3, £3,3 = 3. The third medium is rotated by the angle 6*3 = 1 
with respect to the first one. The angle for the wavevector is chosen as p — 2. The value of the imaginary frequency 
is ^ = IQc. The comparison of the free energy contribution of the retarded calculation (denoted as InCs) and our 
results are given in figure [1] The results again match exactly for this set of parameters as well as for any other one 
might choose. 

This concludes our numerical proof that all the limits with which one can meaningfully compare our result are 
reached exactly within the numerical computation. Again, we have not been able to establish these correspondences 
analytically. The exact expressions are unfortunately extremely convoluted and complicated, thus precluding any 
meaningful and straightforward simplification. 
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FIG. 1: Anisotropic free energy contributions as a function of separation L for the nonretarded case (line) compared to our 
approach (dots). Both free energy expressions coincide numerically to any desired accuracy. See text for details. 
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FIG. 2: Anisotropic free energy contributions as a function of separation L for the retarded case (line) compared to our approach 
(dots). See text for details. 



V. SMOOTH ANISOTROPIC PROFILES 



There is one obviously restrictive condition that is built into the Lifshitz theory of dispersion interactions, which 
underlyes also the approach presented in this work: this is the assumed steplike change in the dielectric permeability 
at the interfaces of bodies that interact across spatially homogeneous media. Several extensions of the theory have 
relaxed the condition of steplike interfaces as well as formulated the dispersion interactions in a non-uniform dielectric 
media [T2| where the response tensors are assumed to vary smoothly and continuously in the z direction. Since 
this extension of the dispersion interactions theory was built on the reformulation of the Lifshitz theory based on 
the algebra of 2x2 matrices, we assume that the case of continuously varying anisotropic response tensors could be 
obtained by a proper generalization of the method derived in this work. 

Indeed it turns out that our method can be applied to the case of smoothly varying response tensors as functions 
of the z variable. As a model, we take an anisotropic dielectric material with the zero frequency tensor components 
ei_i = 5, £2,2 = 10, 63,3 = 20 in its diagonal reference frame. In the coordinate frame where the z direction is 
perpendicular to the translational plane of symmetry, the material is taken to be rotated first by an angle of 1 around 
the X axis, and then by an angle of 2 around the y axis. The corresponding rotated dielectric response tensor is 
denoted by e^. Then, the tensor e^{6) is obtained by further rotating the material by an additional angle 5 around 
the z axis. The model for the dielectric response is then 

( z<0 ; ei 

e{z) = I < z < zo ; cos2 (^) + sin^ (g) e^{S) (65) 

This profile can, for example, be considered as a model for a grain boundary between two grains of an anisotropic 
material p^ . 

We model the continuous dielectric material as a composition of Ni = 100 thinly stratified discrete layers of thickness 
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FIG. 3: Anisotropic free energy density as a function of the rotation angle 5. Solid line denotes the infinite separation free 
energy normalization, the dashed line represents the slowly varying dielectric response normalization. See text for details. 



zq /Ni , where each layer is assigned a constant value of the dielectric response given by equation ()65|) with z taken 
at the middle of the layer. We then apply our general theory to write down the dispersion interaction free energy 
between the two outermost semi-infinite dielectric regions with thinly stratified region in between. The discretization 
can be done with any value of Ni and in the limit approaches the continuum result as was proven exactly for the 



non-retarded isotropic case [14 1 



In order to compute the free energy an assumption needs to be made on both the frequency ^ as well as 
wavevector Q dependence of the dielectric response. For clarity of the model, we assume that the response is constant 
in both ^ and Q until their limiting values of and Qo, at which point the dielectric response becomes 1 throughout 
the material. In the calculations, we choose Qq = 1/zo and = c/zq. 

As has already been mentioned, the free energy calculations require a proper reference point due to the arbitrariness 
of the eigenvector normalizations. Two reference points were taken for calculations at each rotation angle 5. The first 
one was the large empty space separation of the outermost layers as summed up by equation (|43p , and the other was 
by taking a very slowly varying profile of the same shape as given in (j65p but with zq — 30 instead of 1, with the 
relatively low number chosen for numerical stability reasons. The free energy density as a function of the rotation 
angle 5 is shown for both normalization approaches. Fig. [31 The small residual difference between the two methods 
of eigenvector normalization can be attributed to the non-infinite separations in the slowly varying case. 

Our general method for evaluation of dispersion interactions in complicated one-dimensional geometries involving 
non-isotropic dielectric materials obviously predicts the existence of torques, i.e. derivatives of the free energy with 
respect to (5, between the two outermost homogeneous semi-infinite materials. It would thus certainly add to the 
overall energy balance equation of the formation of a grain boundary between two grains of anisotropic crystals, a 
research direction we will not pursue in this paper. 



VI. CONCLUSIONS 



We devised a method that allows us to numerically calculate dispersion interactions in a system composed of an 
arbitrary number of anisotropic and optically active finite dielectric slabs. The result is expressed in the formalism of 
4x4 matrices which, for the most general case, requires a numerical approach. We showed that the results analytically 
reduce to known formulations in the general isotropic case, and that, for a single slab, it gives the same numerical 
result in the nonretarded and uniaxial limits. One of the main strengths of the method is that it is also capable of 
handling smoothly varying media, where the medium is represented as a discretized set of small thickness slabs. 
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